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Abstract 

Several applications exist in which lattice Boltzmann methods (LBM) are used to 
compute stationary states of fluid motions, particularly those driven or modulated 
by external forces. Standard LBM, being explicit time-marching in nature, requires 
a long time to attain steady state convergence, particularly at low Mach numbers 
due to the disparity in characteristic speeds of propagation of different quantities. 
In this paper, we present a preconditioned generalized lattice Boltzmann equation 
(GLBE) with forcing term to accelerate steady state convergence to flows driven by 
external forces. The use of multiple relaxation times in the GLBE allows enhance- 
ment of the numerical stability. Particular focus is given in preconditioning external 
forces, which can be spatially and temporally dependent. In particular, correct forms 
of moment-projections of source/forcing terms are derived such that they recover 
preconditioned Navier-Stokes equations with non-uniform external forces. As an il- 
lustration, we solve an extended system with a preconditioned lattice kinetic equa- 
tion for magnetic induction field at low magnetic Prandtl numbers, which imposes 
Lorentz forces on the flow of conducting fluids. Computational studies, particularly 
in three-dimensions, for canonical problems show that the number of time steps 
needed to reach steady state is reduced by orders of magnitude with precondition- 
ing. In addition, the preconditioning approach resulted in significantly improved 
stability characteristics when compared with the corresponding single relaxation 
time formulation. 
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1 Introduction 

In recent years, the lattice Boltzmann method (LBM) has emerged as an al- 
ternative and accurate approach for computational physics, and, in particular, 
for computational fluid dynamics (CFD) problems [1,2]. It is generally based 
on minimal discrete kinetic models whose emergent behavior, under appropri- 
ate constraints, corresponds to the dynamical equations of fluid flows or other 
physical systems. It involves the solution of the lattice-Boltzmann equation 
(LBE) that represents the evolution of the distribution of particle populations 
due to their collisions and advection on a lattice. When the lattice, which 
represents the discrete directions for propagation of particle populations, sat- 
isfies sufficient rotational symmetries, the LBE recovers the weakly compress- 
ible Navier-Stokes equations (NSE) in the continuum limit. The LBE can be 
constructed to simulate complex flows by incorporating additional physical 
models [3,4]. 

Though it evolved as a computationally efficient form of lattice gas cellular 
automata [5], it was well established about a decade ago that the LBE is actu- 
ally a much simplified form of the continuous Boltzmann equation [6,7]. As a 
result, several previous results in discrete kinetic theory could be directly ap- 
plied to the LBE. This led to, for example, improved physical modeling in var- 
ious situations, such as multiphase fiows [8,9] and multicomponent fiows [10], 
and in an asymptotic theory suitable for rigorous numerical analysis [11]. In 
particular, the latter development has made it possible to study consistency, 
convergence and accuracy of the LBE in a manner similar to the classical nu- 
merical methods for continuum based approaches. As a result of features of 
the stream-and-collide procedure of the LBE such as the algorithmic simplic- 
ity, amenability to parallelization with near-linear scalability, and its ability to 
represent complex boundary conditions and incorporate physical models more 
naturally, it has rapidly found a wide range of appUcations [12,13,14,15,16]. 

Several applications exist where steady state solutions to fiuid fiow problems 
are highly desirable. Examples include magnctohydrodynamic flows and mul- 
tiphase porous media, where one is often mainly interested in investigations of 
their steady state characteristics. On the other hand, the standard form of the 
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LBE is hyperbolic in nature and its solution involves explicit marching in time. 
As a result, it necessarily involves evolving through a transient phase before 
reaching a stationary state. Due to the need to march for many number of time 
steps in this transient phase, it incurs significant computational cost. Another 
important related issue is that the LBE actually represents compressible NSE 
valid at low Mach numbers. Ma. Its deviation from the incompressible NSE, 
which we shall call "compressibility deviations", is independent of grid reso- 
lution. When one intends to simulate close to incompressible flow using LBE, 
such deviations (or the Ma) should be made smaller. This is also desirable 
from a computational viewpoint as the stability regime of the LBE generally 
widens at lower Ma. However, as the Ma is lowered, there is a greater dis- 
parity between the propagation speeds of density perturbations, i.e. the speed 
of sound, and the convection speed of the fluid. In a hyperbolic system, the 
numerical domain of influence should encompass the physical domain [17], re- 
quiring resolution of the time scales of the fluid motion. As a result, computing 
lower Ma flows further compounds the issue and requires a larger number of 
time steps to achieve steady state convergence. 

In recent years, several approaches have been proposed to improve the con- 
vergence rate of the LBE to steady state. These include, in one category, 
reformulations of the LBE to time-independent versions that can be solved 
as a linear system [18,19] and a finite-difference time-independent version 
solved by a multigrid method [20]. In another, they involve adding an artifi- 
cial body force to the time-dependent LBE [21], constructing an implicit LBE 
in a finite-difference or finite-element formulation that allows taking larger 
time steps [22,23,24], or by using a non-linear form of multigrid solver with 
a non-linear LBE time stepping scheme [25]. All these schemes can signifi- 
cantly improve convergence rates, but at the cost of increased complexity as 
compared with the standard LBE. 

On the other hand, Guo et al. [26] proposed an alternate approach to reduce 
the number of time steps necessary for steady state convergence by applying 
preconditioning to the LBE, while maintaining its simplicity. The essential 
principle of this approach, which was originally developed for general hyper- 
bohc schemes by Turkel and others, is as follows [27,28,29,30]. At low Ma, 
in explicit formulations, there is a disparity in propagation speeds of density 
perturbation and fiuid convection. This is formally characterized by higher 
values of condition number, which is defined as the ratio of the fastest to the 
slowest speeds of propagation, or equivalently, the ratio of the maximum to 
minimum eigenvalues of the hyperbolic system, and is inversely proportional 
to Ma. By applying a preconditioner, the speeds of propagation of various 
quantities can be made closer to one another. This can be achieved only at 
the cost of sacrificing the temporal accuracy of the solutions, which in any 
case is not very important as the chief interest is in obtaining steady state 
flow characteristics. Guo et al [26] achieved this in the context of LBE by 



3 



applying a preconditioning parameter that modifies the equihbrium distribu- 
tion function in its collision model. Its emergent behavior is a preconditioned 
compressible NSE with reduced stiffness and hence significantly reduces the 
number of time steps to reach steady state. 

All the preconditioning approaches for LBM mentioned above employ the 
single relaxation time (SRT) model [31] to represent the effect of particle col- 
lisions, with the exception of a recent work that adopts a different approach to 
preconditioning a general form of the LBE than considered here [32]. A com- 
monly used form, the SRT-LBE involves relaxation of particle distributions 
to their local equilibria at a rate determined by a single parameter [33,34]. 
On the other hand, an equivalent representation of distribution functions is in 
terms of their moments, such as various hydrodynamic fields including den- 
sity, mass flux, and stress tensor. The relaxation process due to collisions can 
more naturally be described in terms of a space spanned by such moments, 
which can in general relax at different rates. This forms the basis of the gener- 
alized lattice-Boltzmann equation (GLBE) based on multiple relaxation times 
(MRT) [35,36,37]. By carefully separating the time scales of various hydro- 
dynamic and kinetic modes through a linear stability analysis, the numeri- 
cal stability of the GLBE or MRT-LBE can be significantly improved when 
compared with the SRT-LBE, particularly for more demanding problems at 
high Reynolds numbers [36]. The MRT-LBE has also been extended for mul- 
tiphase flows [38,39,40,41,42], and, more recently, for magnet ohydro dynamic 
problems [43], with superior stability characteristics. It has also been used 
for LES of a class of turbulent flows [44,45,46]. It is known that for a given 
grid resolution and Reynolds number, the standard LBM based on the SRT 
model becomes less stable as Ma is lowered due to the relaxation time becom- 
ing smaller [47,26]. Since the preconditioning is mainly intended to accelerate 
steady state convergence at lower Ma, it is also important to stabiUze the 
computations, which can be optimally achieved by using the MRT-LBE. 

Another consideration is how to precondition the LBE in the presence of ex- 
ternal forces. While Guo et al. [26] suggest a way to precondition a particular 
form of forcing term, details on preconditioning general forms of spatially 
and temporally varying forcing terms are lacking. Such forms are important 
in many situations including magnetohydrodynamic (MHD) flows, where the 
Lorentz force impressed on the fluid can vary spatially and temporally, and 
multiphase flows represented by mean-field models, and buoyant fiows. More- 
over, previous studies were limited to a narrower class of two-dimensional (2D) 
flows, largely in the absence of any body force. In addition, dynamics of flow 
of complex fluids is generally represented by a system of LBE, typically with 
one LBE representing the fiow fields and another one characterizing the evo- 
lution of other physical processes occurring within the fiuid. For example, for 
MHD fiow, we have one LBE to represent the fiuid fiow and another one for 
the magnetic induction equation. Similarly, in the case of multiphase flows. 



4 



we have two sets of LBE - one for the fluid dynamics and another one for the 
dynamics of an order parameter that distinguishes the phases. Thus, it is also 
important to extend the preconditioning to such systems of LBE. 



It is important to note that preconditioning a system of LBE formally improves 
the condition number of its equivalent macroscopic system. For example, in 
the context of MHD, preconditioning a system of LBE actually improves the 
condition number of the equivalent system consisting of the NSE and the 
magnetic induction equation. As a result, while the convergence rate of the 
LBE scheme, which is typically associated with an exponent, is unchanged, 
the prefactor of the convergence rate is modified by preconditioning. In effect, 
the number of time steps needed to reach a steady state representation of the 
equivalent macroscopic system is significantly reduced. 



The primary objective of this paper is then to develop a preconditioning 
method for the MRT-LBE with general forms of forcing terms representing 
non-uniform forces to accelerate convergence to steady state flows. In this re- 
gard, we derive expressions for preconditioned equilibrium moments that gives 
rise to the linear viscous and non-linear convective behavior of a fluid. Based 
on a Chapman- Enskog multiscale analysis [48], we also derive correct func- 
tional forms of the moment projections of source/forcing terms corresponding 
to spatially and temporally dependent variation of forces, which avoids dis- 
crete lattice artifacts. A limiting case of the source terms for the SRT-LBE 
will also be presented. To illustrate the use of preconditioning for a system 
of LBEs, we derive a preconditioned lattice kinetic model for MHD, and also 
provide a simple approach to attain lower values of magnetic Prandtl num- 
ber at steady state, which is important for simulating liquid metal flows. We 
illustrate the advantages of these approaches for a set of canonical problems, 
particularly in three- dimensions (3D). In doing so, we also present some new 
results with shear driven MHD flows. It may be noted that the approach pre- 
sented here, though illustrated for MHD problems, may be readily extended 
to develop preconditioning to a system of MRT-LBEs for a variety of other 
problems. 



This paper is organized as follows. After a brief description of the generalized 
lattice- Boltzmann equation with forcing term in Sec. 2, in Sec. 3 we present 
a derivation of the preconditioned GLBE with forcing term in both 2D and 
3D. The corresponding preconditioned form of lattice kinetic equation for 
magnetic induction is discussed in Sec. 4. Some canonical examples simulated 
using preconditioned LBM are discussed in Sec. 5. Finally, the summary and 
conclusions of this paper are provided in Sec. 6. 
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2 Generalized Lattice Boltzmann Equation with Forcing Term 



The lattice-Boltzmann method computes the evolution of distribution func- 
tions as they move and collide on a lattice grid. The collision process considers 
their relaxation to their local equilibrium values, and the streaming process de- 
scribes their movement along the characteristics directions given by a discrete 
particle velocity space represented by a lattice. Typical lattice velocity models 
include the two-dimensional, nine velocity (D2Q9) and the three-dimensional, 
nineteen velocity (D3Q19) models [33], which are considered in this paper. 
The particle velocity corresponding to the D2Q9 model may be written as: 



= < 



(0,0) a = 

(±1,0),(0,±1) a = l,---,4 
(±1,±1) q; = 5,---,8 



and for the D3Q19 model: 



(0,0,0) q; = 

(±1,0,0),(0,±1,0),(0,0,±1) a = l,---,6 
(±1, ±1, 0), (±1, 0, ±1), (0, ±1, ±1) a = 7, • • • , 18. 



(2) 



The GLBE computes collisions in moment space, while the streaming process 
is performed in the usual particle velocity space [37] . The form of the GLBE 
considered here [46] also computes the forcing term, which represents the effect 
of external forces as a second-order accurate time-discretization, in moment 
space [39,46]. Wc use the following notation in our description of the procedure 
below: In particle velocity space, the local distribution function f, its local 
equilibrium distribution f^'^, and the source terms due to external forces S 
may be written as the following column vectors: f — [fo, fi, f2, ■ ■ ■ , fb]\ = 
[fo", ft", /I', ■ ■ ■ , and S = [^0, S,, S^, . . . , S,]\ where b is the number of 
non-zero discrete velocity directions for a given lattice model. Thus, 6 = 8 
and 6 = 18 for D2Q9 and D3Q19 models, respectively. Here, the superscript 
t represents the transpose operator. 

In particular, the form of the source terms in particle velocity space are ob- 
tained from the expression used in the discrete velocity Boltzmann equation 
— ^ I p ■ fa by approximating it to —'Pip ■ ^ fn'^'^ [49] and fur- 

ther simplifying by neglecting terms of the order of 0{Me?) or higher to 

get [39] Sa = Wa [3 -Tt) + 9{e^-lt)-e'a] -P where Z^'^'*^ = ^^{1 + 
3e^ • it + 9/2(e^ • if )^ — 1/2 if ^} is the local discrete MaxweUian truncated to 
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O(Ma^) [33]. Here, Wa is a weighting factor [33], p and it are the local fluid 
density and velocity, respectively, and 'P is the external force field, whose 
Cartesian components are F^, Fy and F^. 

The moments f are related to the distribution function f through the relation 
f = Tf where T is the transformation matrix. Here, and in the following, 
the "hat" represents the moment space. The transformation matrix T is con- 
structed such that the collision matrix in moment space A is a diagonal matrix 
through A = TAT^^, where A is the collision matrix in particle velocity space. 
The elements of T are obtained in a suitable orthogonal basis as combinations 
of monomials of the Cartesian components of the particle velocity through 
the standard Gram-Schmidt procedure, which are provided by Lallemand and 
Luo [36] and d'Humieres et al. [37] for 2D and 3D lattice models, respectively. 
Similarly, the equilibrium moments and the source terms in moment space 
may be obtained through the transformation f^^ = Tf^*, S = TS. The com- 
ponents of moment-projections of these quantities are: f = fo, fi, f2, ■ ■ ■ , fb , 
feq ^ J-eq^ Jeq ^ jeq ^ • • • , /ft^^ \ and S = Sq, Si, §2, ■ ■ ■ , Sb ^ The exprcssious 
for these quantities are provided in Appendix A for both D2Q9 and D3Q19 
models. 

The solution of the GLBE with forcing term can be written in terms of the 
following "effective" collision and streaming steps, respectively: 

f(^,i) =f(^,i)+'CC7(^,i), (3) 



and 



(4) 



where the distribution function f = {fa}a=o,i,...,b is updated due to "effective" 
collisions resulting in the post-collision distribution function f = {fa}a=o,i,...,b 
before being shifted along the characteristic directions during the streaming 
step. The change in distribution function due to collisions as a relaxation 
process and external forces is represented by -cc^, and following Premnath et 
al. [46] it can written as 



-A (f - f^«) + (l--A 



(5) 



where X is the identity matrix and A = diag(so5 Si, . . . , S;,) is the diagonal 
collision matrix in moment space. Also, here and henceforth, f = f("af,t), 

feq = feg(^^ ^) ^ud S = S(^, t). 
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It may be noted that Eqs. (3) and (4) are obtained from the second-order 
trapezoidal discretization of the source term in the GLBE [39], viz., fa{^ + 
e^St,t + St) - = -J2pKf3 [fpi'^,t) - /|'^(x',t)] + ifa where = 

1/2 af , t) + Sa{'^ + e^St, t + dt)] St, which is made effectively time-explicit 
through a transformation f ^ = fa — ^/'^SaSt [49], and then dropping the 
"overbar" in subsequent representations for convenience. Subsequently, both 
the collision and source terms are represented in the natural moment space of 
GLBE. The second-order discretization provides a more accurate treatment 
of source terms, particulary in correctly recovering general forms of external 
forces in the continuum limit without spurious terms due to discrete lattice 
effects [50], and its time-explicit representation faciliates numerical solution 
in a manner analogous to the standard LEE. 

Now, some of the relaxation times in the collision matrix, i.e. those corre- 
sponding to hydrodynamic modes can be related to the transport coefficients, 
such as the bulk and shear viscosities. The rest of the relaxation parameters, 
i.e. for the kinetic modes, can be chosen through a von Neumann stabihty 
analysis of the hnearized GLBE [36,37]. See also Appendix A for more details. 

Once the distribution function is known, the hydrodynamic fields, i.e., the 
density p, velocity Tt, and pressure p can be obtained as follows: 



b ^1 
P^^fa, y = P^ ^Yl + 7:^St, p = cIp, (6) 

a=0 a=0 ^ 



where, Cg = c/ \/3 with c = 6x/ St being the particle speed, and 6x and 6t are 
the lattice spacing and time step, respectively. 

The computational procedure for the solution of the GLBE with forcing term 
is optimized by fully exploiting the special properties of the transformation 

matrix T: these include its orthogonality, entries with many zero elements, 
and entries with many common elements that are integers, which arc used 
to form the most common sub-expressions for transformation between spaces 
in avoiding direct matrix multiplications [37]. For details, we refer the reader 
to Ref. [46]. As a result of such optimizations, the additional computational 
overhead when GLBE is used in lieu of the popular SRT-LBE is small, typically 
between 15% — 30%, but with much enhanced numerical stability that allows 
maintaining solution fidelity on coarser grids and also in simulating flows at 
higher Reynolds numbers. 
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3 Preconditioned Generalized Lattice Boltzmann Equation with 
Forcing Term for Fluid Flow 

As noted earlier, computation of flows at low Ma using the standard LBE can 
be slow to converge to steady state due to the condition number of its equiv- 
alent NSE being large, as it is inversely proportional to Ma. Moreover, for a 
given Reynolds number, there is a limit on how low Ma can be before numerical 
stability problems result, as the relaxation time in the standard LBE, r, can 
become very close to 0.5 when Ma is made smaller. Preconditioning effectively 
reduces the disparity in propagation speeds of density perturbation and fluid 
convection, or improves the condition number of the equivalent NSE being 
simulated. The use of the GLEE or MRT-LBE improves numerical stability 
by appropriately tuning the relaxation times of the non-hydro dynamic kinetic 
or ghost modes through a von Neumann stability analysis. We now present 
the preconditioned generalized lattice Boltzmann equation with forcing term. 

Several factors need to be considered in preconditioning the GLBE. The 
streaming step in the GLBE is a Lagrangian free-flight or propagation process 
from one lattice node to another node. The collision process is a relaxation 
step that contains linear, faster density propagation process and slower viscous 
momentum transfer process, and non-linear fluid convcctivc process. They are 
individually described in moment space and their separate effects or contribu- 
tions need to be properly preconditioned. Also, careful consideration should 
be given to the preconditioning of the forcing terms in moment space, as their 
contributions, depending on the moment, vary widely, from simple Cartesian 
component of external forces to work due to such forces. In particular, as 
noted in Appendix A, the moment projections of forcing terms are functions 
of external force flelds and velocities, and their products. Hence, care needs to 
be exercised in properly preconditioning individual components of the forcing 
terms corresponding to hydrodynamic and kinetic or ghost modes. As in Guo 
et al. [26], we introduce a preconditioning parameter 7, with < 7 < 1, on the 
GLBE with forcing term. It may be noted that setting 7 equal to 1 reduces 
to the standard form without preconditioning, while 7 < 1 improves the con- 
dition number of the equivalent NSE system of the GLBE. By performing a 
Chapman- Enskog analysis on such GLBE, its preconditioning can be properly 
constructed such that it recovers the corresponding preconditioned compress- 
ible NSE in the continuum limit. The details of this procedure carried out for 
the D2Q9 model is presented in Appendix B, which can be extended to other 
lattice models. 

The preconditioned GLBE with forcing term can be written in terms of the 
following "effective" collision and streaming steps, respectively: 




(7) 
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and 

faC^ + ~^Jt,t + 5t) = fa{'t,t), (8) 

where zu* represents the change in distribution function due to preconditioned 
coUisional relaxation and forcing terms due to external forces. It can be written 
as 



Tu*{-^, t) = T-^ -A* (f - f"«'*) + (2^ - ^A*) S* 



(9) 



Here, I is the identity matrix, A* is the preconditioned diagonal collision 
matrix in moment space, f*^^'* is the preconditioned equilibrium moments and 
S* is the preconditioned moment projections of source terms due to external 
forces. Here, and in the following, the superscript "*" denotes preconditioned 
variables. 

The preconditioning of the components of the equilibrium moments, 

f^^'* = [/o''*,/r*,/i''*,---,/r'f (10) 

which are functions of the conserved moments, can be performed by analyzing 
the GLBE using the Chapman- Enksog expansion, as in Appendix B. The 
components of F'^ * can be written for the D2Q9 model as: /q*'* = p, /f*'* = 

_A feq,* _ ■ req,* _ eg,* — _A feq,* _ eq,* _ Ux-jy) feq,* _ eq,* _ j^jy 
Jx-)J5 JyiJe ~% Jy>J7 — Pxx 'yp 'J8 ~ fxy ■ 

The definition of the components of the equilibrium moments are provided in 
Appendix A. 

For the D3Q19 model, they become: /q'^'* = p, /f'^'* = e^''-* = -llp+ 19-^, 

feq,* _ 2,eq,* _ q„ _ 11 3 ■ j feq,* _ ■ feq,* _ eq,* _ _2 • feq,* _ 
J2 — c —op ^^)J3 — JxtJa — Hx — 3Jx:J5 — 

■ feq,* _ eq,* _ _2 • feq,* _ ■ feq,* eq,* _ _2 • feq,* _ q„eg,* _ 

3£- j ■ j 



IP 



feq,* _ n^eq,* _ q /_l„eg,*\ feq,* _ eq,* _ [^v ^i] feq,* _ , 



eq,* 

WW 



_1 eg,* ieq,* _ eq,* _ Jxjy feq,* _ eg,* _ JyJz ieg,* _ eg,* _ jxjz feq,* _ 
2"uiui'^13 — Pxy iJ 14 — Pyz ^p i J 15 — Pxz jp ' J 16 

0, /rr = 0, iri'* = o. 

A general observation is that only the non-linear terms in the components 
of the equilibrium moments are preconditioned by the parameter 7. This is 
consistent with the argument that the hydrodynamic convective effects, which 
are non-linear, emerge from relaxation process during collisions should be con- 
tained in these terms; they should be preconditioned to match the faster prop- 
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agation of density pcrtTirbations, which arc represented by linear terms in the 
equihbrium moments. It may be noted that an alternative approach to pre- 
conditioning the equilibria has been proposed recently [32]. 

The preconditioned components of the source terms 



S* 



Q* C* C* 

Oq, O^, O2, . . . , 



fill 



can be written, for the D2Q9 model as: = 0,SI = y^y) , S* = 

(FxUx+FyUy) O* Fx Q* Fx Q* Fy O* Fy Q* Q {FxUx—FyUy) Q* 

~0 ^3 ) '-'3 — ~) '-'4 — ~~i>~'5 — ~i'~'6 — ~~i'~'7 — ^ ' '-'8 — 



The corresponding components of S* for the D3Q19 model are: = 0, = 

OQ (FxUx+FyUy+FzUz) Q* 1 {FxUx+ FyUy + F^U z) Fx Q* 2Fx Q* Fy 

■JO -y2 J '-'2 — -,2 • — J : '^4 — 3 7 ' '-'5 ~ -y ' ^6 ^ 
2Fy qst: Fz c* _2_F^ c* r) {2K, n J tly-F.z'Lz) 

'-'10 ~ 72 ; ^11 — ^ ^2 , 

Q* _ (FyUy-FzUz) Q* _ (F^Uy+FyUx) Q* _ (FyUz+FzUy) Q* _ (F^Uz+FzUx) Q* 

■ y2 T '~>13 — 72 T '-'U — 72 ) '-'15 ~ 72 )*-'l6~ 

The preconditioning of the moment projections of the source terms may also 
be compactly written as 

S* = VsS, (12) 



where 



^ 1 1 1 1 1 1 1 M 

7^5 = diag 1, — (13) 

V 0^ 0^ 0^ 0^ 0^ 0^ 0^ 0^ / 



for the D2Q9 model, and 



_ 1 11111111 1 1 1 1 1 1 .W.^N 

' <S Q-iag I i, r-, , , , , , , r, r-, r, i, i, i I (^-LT:j 

V 0' 0' ^ ^ ^ ^ ^ ^ J 



for the D3Q19 model, where the components of the unpreconditioned source 
terms S are given in Appendix A. 

Clearly, the external forces have a first-order effect on the convective motion 
of the fluid, and thus to "condition" the moments linearly influenced by such 
forces, the moment projections need to be preconditioned by the inverse of 7. 
On the other hand, other moments are effected by the external forces at second 
order. These include the "work" contribution by their interaction with the fluid 
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motion on the moment corresponding to kinetic energy (see Appendix A). 
Similarly, the moment projections of the source terms for the momentum flux 
tensors have second-order influence. In general, the Chapman-Enskog analysis 
reveals that all higher order moments that involve non-linear effects from 
interaction of external forces and fluid motion are much slower than the fluid 
motion itself and needs to be preconditioned by the inverse of the square of 
the preconditioning parameter, i.e. 1/7^ (see Appendix B). 

For the preconditioned collision relaxation time matrix, 

A* = diag(s*,st,...,s*), (15) 



some of the relaxation times s*, i.e. those corresponding to hydrodynamic 
modes, can be related to the transport coefficients. The rest, i.e. those for the 
kinetic modes, can be chosen through a von Neumann stability analysis of the 
hnearized GLBE [36,37]. For the D2Q9 model, we have i = 3^ + i where 

C is the bulk viscosity, and Sy = Sg = s*, where ^ = 3^ + |, with u being 
shear viscosity. For the kinetic modes, we have [36] si = 1.63, ^2 = 1.14 and 
= si = 1.92. On the other hand, for the D3Q19 model [37], we have, for 

the hydrodynamic modes, = | - + ^, Sg = s]']^ = 5*3 = 5*4 = s^g = s*, where 

= 3^ + 1 and for the kinetic modes [37], = 1.19, S2 = •^lo = •^12 = 1-4, 
sl = sl = sl = 1.2 and s^g = s^y = s^g = 1.98. 

The hydrodynamic fields, i.e., the density p, velocity Tt and pressure p ob- 
tained from the solution of preconditioned GLBE, satisfy the equivalent pre- 
conditioned compressible NSE (see Appendix B), and can be written as 

P^Y^fa, J = P^l? = I^/a"e'a + -— 5t, P=cfp, (16) 
a=0 a=0 ^ 7 



where, c* = ^/7Cs with = c/ \/3- The preconditioning of the GLBE effec- 
tively reduces the speed of sound by a factor ^/j. As a result, the disparity be- 
tween the propagation speed of density perturbation and that of fiuid motion 
is decreased by decreasing the parameter 7. Moreover, the "effective" Mach 
number after preconditioning is Ma* = u/c* = Ma/^{j). It may be noted 
that a Chapman-Enskog analysis of the GLBE carried out in Appendix B) 
also shows how the evolution of kinetic modes, in addition to the hydrody- 
namic modes of interest, are affected by preconditioning. It is evident that the 
structure of the preconditioned GLBE with forcing term is very similar to that 
without preconditioning, involving only local scaling of the equilibrium mo- 
ments, the moment projections of source terms and the relaxation matrix. As 
a result, the optimized computational procedure for GLBE with forcing term 
described in the previous section can be fully exploited for the preconditioned 
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version. 



3. 1 Limiting Form: Preconditioned SRT-LBE with Forcing Term 



When all the relaxation parameters are set to the same constant, i.e. = 
1/r*, we arrive at the SRT-LBE, which can be conveniently written as the 
following collision and streaming steps, respectively, where both steps are ex- 
pressed in particle velocity space: 

t) = {fa - rjn + (i - ^) s:st, (i?) 

where = fa{x',t), f^"'* = and S*^ = S*{lt,t), and 

+ e^5t, t + 5t)^ fa(-^, t). (18) 



Here, /^^'* is the preconditioned equilibrium distribution 



feq,* 
J a 



27^ " ' 27 



(19) 



The hydrodynamic fields can be obtained from the distribution functions in 

the same manner as before, i.e. from Eq. (16). One important consideration 
is in obtaining the correct expression for the corresponding preconditioned 
source terms. They can be obtained simply by an inverse transformation of 
the moment projections of preconditioned source terms from Eq. (12). That 
is, 

S* = T-^S* = T-^VsS. (20) 



Explicit evaluation of this equation, Eq. (20), yields 



7 



+ 9- 



3q • li ) e a ' 



(21) 



which is the desired expression for the source term of the SRT-LBE with 
preconditioning. It may be noted that it is essential to maintain the above 
form of preconditioned source term to correctly recover the corresponding 
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preconditioned hydrodynamic behavior. On the other hand, for example, if, 
one naively sets 

S: = ^ ^ i, (22) 

a Chapman-Enskog analysis of the resulting SRT-LBE (not exphcitly shown 
here, for brevity) yields macrodynamical equations with non-vanishing spuri- 
ous terms (with 7 < 1). The i-th Cartesian component of these extra spurious 
terms to the corresponding preconditioned momentum equations turns out to 

be 

Extra Termj = dj 

These terms can indeed dominate with particularly strong preconditioning at 
lower 7, when (7 — l)/7^ can become very large, particularly for spatially and 
temporally dependent external forces. For example, simulation of MHD prob- 
lems, where Lorentz forces can vary both in space and time, using a precon- 
ditioned SRT-LBE with Eq. (21) yielded accurate results, but with Eq. (22), 
it resulted in grossly wrong behavior. This stresses the critical importance of 
properly preconditioning forcing terms, as the temporal change in the effect 
of the external forces on various physical processes during coUisional relax- 
ation are different. In this regard, analysis of their contributions in moment 
space, as shown above, is particularly revealing: the individual contributions 
of the external forces spanned in the moment space need to be separately 
preconditioned depending on the nature of their effects on the moments. 



4 Preconditioned Vector Lattice Kinetic Equation for Magnetic 
Induction 

As an illustration of preconditioning an extended system of LBE for complex 
fluid flows subjected to external forces, we will now discuss preconditioning 
lattice kinetic equations for the magnetic induction equation required for sim- 
ulation of MHD flows. Dellar [51] concluded that a vector formulation of the 
kinetic equation is necessary to properly recover the magnetic induction equa- 
tion and constructed a 2D model to accomplish this, which was extended to 
3D by Breyannis and Valeougeorgis [52] . The GLBE with forcing term is used 
in conjunction with such a lattice kinetic equation, with the latter providing 
the Lorentz force field to the former. 

In addition to the propagation of the density perturbation as sound with speed 



(7-1) 
7^ 



(23) 
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Cs, MHD flows are characterized by the propagation of perturbation of mag- 
netic induction, the so-called Alfven waves. If Bi is the Cartesian component 
of magnetic induction, we can obtain the corresponding Alfven velocity as 
Va^i — Bi/y/pjl, where p and /j, are the density and magnetic permeability, 
respectively. Thus, we can define a local Alfven number Al — VA,i/cs, and 
Dcllar [51] constructed a lattice kinetic equation that recovers the magnetic 
induction equation applicable at low Al, with deviations 0{Af). In this scal- 
ing, 0{VA,i) ~ 0{ui), or 0{A1) fti 0(Ma). Thus, in MHD flows, there is an 
additional disparity between the speed of perturbation of the magnetic induc- 
tion field and the the speed of sound. The condition number, in this case, is 
inversely proportional to the Alfven number, i.e. 0(1/Ai). Preconditioning the 
lattice kinetic equation accelerates its steady state convergence by reducing 
the disparity between such characteristic speeds in MHD flows. 

We now develop a preconditioning formulation for a unified vector lattice ki- 
netic equation for magnetic induction applicable in both 2D and 3D. Unlike 
the case of fiuid flow, which has fourth-order isotropy requirements on lat- 
tice velocity models to correctly recover viscous stress tensor, the magnetic 
induction imposes lower order symmetry requirements. Thus, we need only a 
smaller number of particle velocity directions for magnetic induction e^"*, and 
following previous work, we consider a D2Q5 model 



(0,0) a^O ^^^^ 

(±l,0),(0,±l)a = l,---,4 



and a D3Q7 model 

(0,0,0) a = 



^ a 



(25) 

±1, 0, 0), (0, ±1, 0), (0, 0, ±1) a = 1, • • • , 6 



in 2D and 3D, respectively. The preconditioned lattice kinetic equation can, 
then, be written in terms of the following collision and streaming steps, re- 
spectively: 



9aj 



i^,t)^-j^{g^,-g:r) (26) 



where g^j = gajC^^t) and g^Y = Qaf^'^^'t), and 

gaj{-t + "e'^/t, t + 5t)^ gajC^, t), (27) 

where gaj is the vector distribution function in index notation o; = 0, 1, • • • , fo^. 
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Here, 6^ = 4 and 6^ = 6 in 2D and 3D, respectively. The subscript Roman 
indices i, j, etc., represent Cartesian components of the coordinate directions. 
Assuming the usual summation convention of repeated indices, the Cartesian 
component of the preconditioned vector equilibrium distribution function g^j* 
is given as 




(28) 



where 



W^a=^ (29) 
[l/{2Nm) a = !,■■■, bm 



and 6m = with = 3 and N„i = 4 in 2D and 3D, respectively. 

Here, Bj is the Cartesian component of the magnetic induction and 7^ is 
the preconditioning parameter, with < 7^ < 1. Thus, the preconditioning 
is carried out on the non-linear part of the vector equilibrium distribution, 
Eq. (28), which represents the transport of magnetic induction field by the 
fluid motion. The preconditioned relaxation time is related to the magnetic 
diffusivity of the fluid rim, where 77^ = 1/ (/la) with /i and a being the magnetic 
permeability and electrical conductivity, respectively, and is given as 

= + ^. (30) 



Once the vector distribution function is calculated, the components of the 
magnetic induction Bi and the current density Jj can be obtained by taking 
their zeroth and first moments as: 

bm 

B^ = Y.9ai (31) 

a=0 



and 



1 1 bn. 

^ ^ \ ^ / e.q,*\ 

1^ 'T'm^rn a=0 



(32) 



where eijk is the Levi-Civita or the third-order permutation tensor. It can be 
shown that the preconditioned vector lattice kinetic equation can recover the 
corresponding preconditioned lattice induction equation given as follows (see 
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Appendix C): 

1 1 
dtBi H Vj {ujBi - BjUi) = — dj {rjmdjBi) . (33) 

As shown by Dellar [51], the magnetic induction will remain solenoidal, i.e. 
diBi = 0, provided the initial condition on the magnetic induction satisfies 
the divergence free condition. 

The interaction of the magnetic induction and the current density gives rise 
to the Lorentz force on the fiuid fiow. This force can be written as 

'^Lorentz = ~f X ^ (34) 

and enters as ^ = '^Lorentz in the preconditioned GLBE discussed in the 
previous section. 



4.1 Achieving Low Magnetic Reynolds Number or Magnetic Prandtl Number 
at Steady State 



If L, mq, Bq are the characteristic length, velocity and magnetic induction 
scales, respectively, then we can non-dimensionalize various quantities as ^ 
, t ^ {uo/L)t, u ^ u/uq and B ^ B/Bq, and obtain the non-dimensional 
form of magnetic induction equation, Eq. (C.IO) as 

ImdtBi + Wj {ujBi - BjUi) = dj ' (35) 

where Re^ = UQL/r]m is the magnetic Reynolds number and the corresponding 
dimensionless current density is 7 = (1/ReJ^ X 

Now, in many practical MHD applications, particularly for liquid metals, Re„i 
is relatively very small, or so is the magnetic Prandtl number Pr^, which is 
given by Pim = vjr] = Rem/ Re: Re^n -C 0(1) and Pr^^ -C 0(1). To achieve 
lower J?e^ or Pr^, we need to make rj^. smaller, which, in turn, from Eq. (30), 
means reducing r^. However, true for a typical lattice based method, as 
approaches 0.5, it can cause numerical instability. This situation can be reme- 
died for the steady state situation by considering the following: at steady 
state, we have Vj {ujBi — BjUi) = dj ( j^g djBij, to which we apply a scaling 

factor X X^j {ujBi — BjUi) = dj ( j^g djBi^. This effectively changes Re^ 
to xR^m at steady state. Thus, we can write Rem,eff = xR^m, resulting in 
Pi'm,eff = xP^m- In the context of preconditioning, a lattice kinetic scheme 
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for this "effective" steady state magnetic induction equation can readily be 
constructed. 

Thus, in dimensional form, we need to construct a scaled lattice kinetic scheme 
with preconditioning for the following macroscopic equation 

X 1 
dtBi H Vj {ujBi - BjUi) = — dj {rjmdjBi) . (36) 



The magnetic field satisfying this equation, Eq. (36), can be obtained by solv- 
ing the above preconditioned lattice kinetic scheme, i.e. Eqs. (26) and (27), 
with its associated auxiliary equations, except for the following changes in the 
computation of vector equilibrium distribution and current density: 

faT = ^« + M - B,uA (37) 



and 

1 ^ 11 



5 Results and Discussion 



We will now present investigations of the preconditioned computational ap- 
proach presented in the previous sections by means of a set of canonical exam- 
ples. Unless otherwise stated, all the results will be expressed in the natural 
lattice units of the method, i.e. we use the lattice spacing 5x as the length 
scale and the particle velocity c as the velocity scale (with 5x/ c used to scale 
the temporal quantities). 

First, we simulate the simple classical problem of flow with a fluid viscosity 
v between parallel plates spaced 2L apart and driven by a pressure gradient 
—dp/dx, i.e. plane Poiseuille flow using the preconditioned GLBE with forcing 
term. We consider the domain to be periodic in the streamwise and spanwise 
directions, and thus the pressure gradient is applied as a body force. No slip 
conditions at the walls are specifled by using the half-way or link bounce 
back scheme [53]. For this setup, if ^ = —{dp/dx)i is the driving force, the 
maximum fluid velocity occurring at the center is u^ax — FL'^ /{2pQi'), where 
po is the nominal fluid density. We set L = 32, u = 0.001 and po = Ij and 
apply a pressure gradient such that the maximum velocity is Umax = 0.00051, 
or Ma = 0.0008333. The Reynolds number based on the above velocity and 
L becomes Re — 32.6. Figure 1 shows the convergence history with different 
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Fig. 1. Convergence histories (in log- log scale) of preconditioned MRT-LBE 
with forcing term for simulation of Poiseuille flow with maximum flow velocity 
'^max = 0.00051 and kinematic viscosity v = 0.001; Re = 32.6. 

values of preconditioning parameter 7 (no preconditioning corresponds to 7 = 
1) for this external force driven flow problem with relatively low Ma. The 
convergence to steady state is measured by the second-norm of the residual 
error of the velocity, | \u{t+ 10) — u(t) 1 12. It can be seen that preconditioning the 
GLBE with forcing term dramatically accelerates the steady state convergence 
for this problem, in particular, by more than two orders of magnitude with 
strong preconditioning carried out using 7 — 0.001. 

Next, we simulate plane Poiseuille flow with relatively higher Ma and Re. As 
before, we set L = 64, but use v = 0.005 and apply a pressure gradient such 
that Ma = 0.00222 and Re = 163.8. The convergence history for this problem 
is presented in Fig. 2. Again, a significant reduction in the number of time 
steps for convergence to steady state is achieved through preconditioning. On 
the other hand, it appears that to maintain numerical stability the minimum 
possible value of the preconditioning parameter 7 needs to be higher at higher 
Ma. The computed velocity profile for problem with preconditioning using 
7 = 0.1 compared with the analytical solution is shown in Fig. 3. Excellent 
agreement is seen. 

When all the other parameters are maintained constant, the deviation of the 
computed solution from the analytical solution is related to the ratio z//7, 
which, in turn, is related to the hydrodynamic relaxation time in the precon- 
ditioned GLBE as ^ = 3- -I- I (see Sec. 3). Figure 4 shows the relative error 
in the computed velocity as a function of u/'j. It is evident that the error. 
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Fig. 2. Convergence histories (in log- log scale) of preconditioned MRT-LBE 
with forcing term for simulation of Poiseuille flow with maximum flow velocity 
Umax = 0.0128 and kinematic viscosity u = 0.005; Re = 163.8. 
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Fig. 3. Comparison of computed velocity profile using the preconditioned MRT-LBE 
with analytical solution for simulation of Poiseuille flow with preconditioning pa- 
rameter 7 = 0.1 and kinematic viscosity v = 0.005; Re = 163.8 
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Fig. 4. Error in the computed velocity at the center by using the preconditioned 
MRT-LBE for simulation of Poiseuille flow as a function of the ratio 1//7. 

which remains relatively small, is linearly proportional to this ratio. Thus, by 
maintaining low values of 1^/7, the error can be correspondingly kept small. 

We now investigate the influence of preconditioning the GLBE on numeri- 
cal stabihty for the case of plane Poiseuille flow considered above. This can 
be done by systematically carrying out simulations at different characteris- 
tic parameters, including 7, u and Ma, and determine the threshold at which 
the computations become unstable, i.e. small variations growing exponentially 
with time. The results can be conveniently expressed in the form of a regime 
map or parameter space that delineates stable and unstable parameter sets. 
Figure 5 shows the stability-instability parameter space determined by the 
maximum flow velocity Umax and the preconditioning parameter 7 for differ- 
ent fluid viscosity u. Arrows normal to the curves pointed upwards indicate 
stable parameter space, while those pointed downwards pertain to unstable 
space. This regime map is particularly revealing. First, for a given fluid viscos- 
ity ly, as the maximum velocity or, equivalently. Ma is reduced, lower values of 
the preconditioning parameter 7 can be used, i.e. the GLBE can be strongly 
preconditioned resulting in greater computational gains while maintaining nu- 
merical stability. The fluid viscosity appears to significantly affect the stability 
parameter space. For a given Umax, the minimum possible value of 7 is higher at 
higher values of u. That is, the extent of preconditioning is greater with lower 
fluid viscosities. Interestingly, Fig. 5 also shows that the delineating curve has 
a linear functional relationship between 7 and Umax at higher u, while at lower 
u, the stability envelope is nearly flat with a constant 7 for a wide variation 
of Umax- Thus, in general, the benefits of preconditioning, while maintaining 
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Fig. 5. Parameter space of the preconditioning parameter 7 for stable/unstable 
computations using the preconditioned MRT-LBE for simulation of Poiseuille flow 
at different maximum possible flow velocities Umax smd kinematic viscosities u of 
the fluid. 

numerical stability, is greater at lower Ma and lower v. It may be noted that 
the use of GLBE in lieu of SRT-LBE in the context of preconditioning results 
in significant enhancement of numerical stability, as will be shown later for 
fully 3D problems characterized by complex fluid motions. 

Let us now consider a problem in which we can investigate preconditioning 
a system of LBE, where the external force impressed on the fluid is a also 
a strong function of the fluid motion itself. This is, we consider the classical 
Hartmann flow, consisting of pressure driven flow between two parallel plates 
in the presence of a magnetic field perpendicular to the walls. In addition 
to the Reynolds number Re, this problem is characterized by the Hartmann 
number Ha, which is given by Ha — BqL^J^, where Bq and a are the applied 
magnetic field strength and fluid's electrical conductivity, respectively, and 
other parameters are as previously defined. 

We now simulate this problem by using a vector lattice kinetic scheme precon- 
ditioned by parameter 7^ for magnetic induction (as in Sec. 4), in conjunction 
with the GLBE with forcing term preconditioned by parameter 7 (as in Sec. 3). 
The Lorentz force arising from the interaction of the magnetic induction and 
velocity field is introduced into consistently preconditioned forcing terms. To 
simulate the fiow of liquid metals at low magnetic Reynolds number or 
magnetic Prandtl number Pr^,, we further apply a scaling factor % to the 
preconditioned lattice kinetic scheme (as in Sec. 4.1). 
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Fig. 6. Comparison of computed velocity profile using the preconditioned MRT-LBE 
with analytical solution for simulation of Hartmann flow with preconditioning pa- 
rameters 7 = 7m = 0.05 and kinematic viscosity v = 0.004; Re = 286, Jfa = 71.6. 

Figure 6 shows the computed steady state velocity profile of Hartmann flow 
with He = 286, Ha = 71.6, Pr^ = 1.0 x 10"^. This is achieved by using 
transport coefficients of u = 0.004, and 1] = 0.004, and L = 64, along with 
preconditioning and scaling factors of 7 = 0.05, 7^ = 0.05 and x = 1.0 x 10^^. 
The corresponding computed steady state magnetic induction profile is pre- 
sented in Fig. 7. The fiattening of the velocity profile observed is characteristic 
of MHD flows due to the Lorentz force, with most of the velocity variation 
being conflned to the region very close to the wall, in the so-called Hartmann 
layer. In general, the thickness of Hartmann layer 6h is inversely proportional 
to the Hartmann number {Sh ~ 1/ Ha) . The analytical solution for this prob- 
lem is provided in standard texts on MHD (e.g., Ref. [54]). It is seen that 
these proflles computed using preconditioned system of LBE are in excellent 
agreement with corresponding analytical solutions. 

Let us now investigate the behavior of steady state convergence of this problem 
by applying various levels of preconditioning. Figure 8 shows the convergence 
history obtained at various values of 7^ at a flxed value of 7, i.e. 7 = 0.05. 
Interestingly, preconditioning only the GLBE with forcing term, but not the 
vector lattice kinetic scheme (i.e. with 7^ = 1.0) results in the slowest conver- 
gence. However, as we increase the extent of preconditioning by lowering the 
values of 7^, the number of time steps to reach steady state is signiflcantly 
reduced. The beneflts of preconditioning are greatest when both the precondi- 
tioning parameters are equal to one another, i.e. 7 = 7m = 0.05. Moreover, it 
is also interesting to observe that if the lattice kinetic scheme is more strongly 
preconditioned than that for the GLBE, i.e. 7m < 7, the approach to steady 
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Fig. 7. Comparison of computed induced magnetic field profile using the precon- 
ditioned MRT-LBE with analytical solution for simulation of Hartmann flow with 
preconditioning parameters 7 = 7m = 0.05 and kinematic viscosity v = 0.004; 
Re = 286, Jfa= 71.6. 



Y = Pr =1.0x10"" 

m 

V = 0.004 
ri = 0.004 




Fig. 8. Convergence histories (in log- log scale) of preconditioned MRT-LBE with 
forcing term in conjunction with preconditioned vector lattice kinetic scheme for 
simulation of Hartmann flow with flow preconditioning parameter 7 = 0.05 for 
various values of induction preconditioner parameter 7m! 

Re = 286, Ha = 71.6. 
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Fig. 9. Convergence histories (in log-log scale) of MRT-LBE with forcing term with- 
out preconditioning in conjunction with preconditioned vector lattice kinetic scheme 
for simulation of Hartmann flow with flow preconditioning parameter 7 = 1.0 for 
various values of induction preconditioner parameter ■7m! 

Re = 286, Ha = 71.6. 

state becomes slower as compared with the case with 7^ = 7. This is consis- 
tent with the scaling 0{AI) ~ 0(Ma), that is, the propagation speeds of both 
magnetic field and velocity field by the fluid convective motion occur at similar 
time scales. Both these are slower than the speed of density perturbation and, 
thus, preconditioning the terms representing these two physical processes by 
the same magnitudes would result in the fastest steady state convergence. 

This notion is further illuminated by considering cases in which the GLBE is 
used without preconditioning (7 = 1), while the lattice kinetic scheme is solved 
with different values of the preconditioning parameter 7^. The results with 
these tests are plotted in Fig. 9. It shows that the convergence rate actually 
becomes slower if only one of the LBEs is preconditioned, no matter what the 
value of the preconditioning parameter is used, when the other LBE is used 
without preconditioning. Thus, preconditioning should be done for both the 
LBE and at the same levels. This is further corroborated by considering a 
series of cases with 7 = 7^, as shown in Fig. 10, with lower values for these 
parameters providing more rapid convergence to steady state. 

It may be noted that, in a recent work, we studied a series of problems with 
very thin Hartmann layers by introducing stretched grids through the Roberts 
boundary layer transformation [55] by means of an interpolation-supplemented 
streaming step in the LBE framework [43]. In particular, using this modified 
LBE, we simulated Hartmann fiow with Ha as high as 10, 000 in very good 
comparison with corresponding asymptotic analytical results, leading to very 
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Fig. 10. Convergence histories (in log-log scale) of preconditioned MRT-LBE with 
forcing term in conjunction with preconditioned vector lattice kinetic scheme for 
simulation of Hartmann flow with equal values of flow (7) and induction (7m) pre- 
conditioning parameters; Re = 286, Ha = 71.6. 



significant reduction in the number of grid nodes as compared to using stan- 
dard LBE using uniform grids [43]. 

We now present some applications of preconditioned LBE for simulation of 
3D wall bounded shear flows of electrically conducting fluids such as liquid 
metals mediated by magnetic fields. In this regard, we consider a cubic cavity 
of side length W containing liquid metal which is driven by its top lid moving 
at velocity Uq. An external magnetic field Bq is applied parallel to the lid 
surface and perpendicular to its direction of motion. A schematic of this fiow 
problem is shown in Fig. 11. 

We consider the top lid moving with a velocity Uq = 0.0235 imparting shear 
on the fiuid of viscosity u — 0.015 contained in a cavity discretized by 64^ 
grid nodes, such that the Reynolds number is 100. Initially, we clarify the 
advantage of preconditioning for this highly 3D fiow even for a simpler case 
that does not involve the application of magnetic field. The velocity boundary 
conditions at the walls, including the top moving lid, were imposed by using 
a link bounce back scheme [53]. For the moving wall, this scheme adds con- 
tributions due to appropriate momentum to the distribution functions. The 
convergence histories in the absence of an external magnetic field for different 
values of 7 are shown in Fig. 12. It is seen that significant reduction in the 
number of time steps to reach steady sate, is obtained with preconditioning 
for this problem. 
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Fig. 11. Schematic illustration of a three-dimensional (3D) cubical cavity with its 
top lid moving at velocity Uq in the presence of an external magnetic field Bq, 
applied parallel to the lid surface and perpendicular to its direction of motion. 




Time step 



Fig. 12. Convergence histories (in log-log scale) of preconditioned MRT-LBE with 
forcing term for simulation of flow in a 3D cubical driven cavity in the absence 
of an applied magnetic field with lid velocity Uq = 0.0235 and kinematic viscosity 
u = 0.015 at different values of flow preconditioning parameters 7; Re =100, grid 
resolution is 64^. 
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Fig. 13. Convergence histories (in log-log scale) of preconditioned MRT-LBE with 
forcing term for simulation of flow in a 3D cubical driven cavity in the presence 
of an applied magnetic field with lid velocity C/q = 0.0235 and kinematic viscosity 
V = 0.015 at different values of preconditioning parameters, with 7 = 7^; fie =100, 
ffa = 45, grid resolution is 64^. 



We now consider the case involving the application of an external magnetic 
field such that the Hartmann number is 45, which is obtained by setting 
r\ = 0.015, and the magnetic Prandtl number is 5.625 x 10~^. We consider 
the induced magnetic fields at far-off distances outside of the cavity to be zero 
as our boundary condition. This is achieved by considering a larger computa- 
tional domain for the magnetic induction field encompassing the cavity walls. 
On these extended computational boundaries, we implement a zero induced 
field condition through an extrapolation method [56] applied to the vector 
distribution functions. The corresponding convergence histories for this 3D 
MHD fiow problem are shown in Fig. 13. Again, a significant reduction in the 
number of time steps to reach steady state is achieved through precondition- 
ing. Moreover, to put things in perspective, when the SRT-LBE was employed 
for the same grid resolution as above, the simulations were found to be stable 
only for u > 0.166. On the other hand, with GLBE, as indicated above, we 
could use a much lower value for 0.015 while maintaining numerical stability. 
Thus, for this problem, by using the preconditioned GLBE rather than the 
preconditioned SRT-LBE, the numerical stability is enhanced by almost an 
order of magnitude. 

We will now investigate the accuracy of preconditioned LBE for this problem. 
Figures 14, 15 and 16 show the computed velocity profiles for the cases with 
Ha — (i.e. no magnetic field) and Ha — 45 and compared with recent results 
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Fig. 14. Comparison of computed velocity profile u along the line x = W/2 and 
y = W/2 using the preconditioned MRT-LBE with finite-difference solution of NSE 
(Morley et al. (2004)) for simulation of 3D cubical lid-driven cavity flow with and 
without magnetic field, i.e., Ha = and Ha = 45, respectively; Re = 100. 
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Fig. 15. Comparison of computed velocity profile u along the line x = W/2 and 
z = W / 2 using the preconditioned MRT-LBE with finite-difference solution of NSE 
(Morley et al. (2004)) for simulation of 3D cubical lid-driven cavity flow with and 
without magnetic field, i.e.. Ha = and Ha = 45, respectively; He = 100. 

from simulations carried out using finite-difference method (FDM) involving 
the solution of the Navier-Stokes equations [57]. The presence of Lorentz 
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Fig. 16. Comparison of computed velocity profile w along the line y = W/2 and 
z = W/2 using the preconditioned MRT-LBE with finite-difference solution of NSE 
(Morley et al. (2004)) for simulation of 3D cubical lid-driven cavity flow with and 
without magnetic field, i.e., Ha = and Ha = 45, respectively; Re = 100. 

forces appears to significantly influence the characteristics of fluid motion in 
this 3D problem. In particular, the velocity profile appears to be markedly fiat- 
tened by the presence of magnetic field. The computed results are in excellent 
agreement with the FDM. It was noticed that the velocity profile bounded 
by the Hartmann walls (i.e. those perpendicular to the direction of Bq, see 
Fig. 11), is somewhat sensitive to the grid resolution employed. Morley et 
al. [57], who studied this problem using different numerical schemes with dif- 
ferent grid resolutions, also observed such effects. In this work, we find that 
by further refining the grid, i.e. by doubling the number of grid nodes in each 
direction, the computed solution converges to the FDM results. 

For the sake of completeness, we will now consider 3D MHD driven cavity fiow 
with the magnetic field applied in a different manner, i.e. Bq perpendicular to 
the lid surface and its direction of motion as shown in Fig. 17. 
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Fig. 17. Schematic illustration of a three-dimensional (3D) cubical cavity with its 
top lid moving at velocity Uo in the presence of an external magnetic field Bq, 
applied perpendicular to both the lid surface and its direction of motion. 
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Fig. 18. Effect of magnetic Prandtl number Pr^ on computed velocity profile u 
along the line x = W/2 and y = W/2 for simulation of 3D cubical lid-driven cavity 
in the presence of magnetic field perpendicular to the lid surface and its direction 
of motion; Re = 100 and Ha = 45. 

We will consider the effect of magnetic Prandtl number for this case by 
employing different values of the scaling factor x in the preconditioned LBE. 
Figures 18, 19 and 20 show the velocity profiles along different directions 
for He = 100, Ha — 45 with 7 = 7^ = 0.1 for two values of Pr^, i.e., 
Pr^ = 5.6 X 10~^ and Pr^ = 1.0, with the former corresponding to liquid 
metal. It is noticed that the values of Pr^ appears to strongly modulate 
the velocity field for those cases in which they are bounded on both sides by 
stationary walls (Figs. 19 and 20). Moreover, the direction of the application 
of magnetic field appears to have a profound infiuence on the fiow field. In 
particular, in contrast to the previous case, we notice the presence of wall jet 
like features when Bq is perpendicular to both the surface of the top lid and 
its direction of motion (see Fig. 21). 
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Fig. 19. Effect of magnetic Prandtl number Pr^ on computed velocity profile u 
along the line x = W/2 and z = W/2 for simulation of 3D cubical lid-driven cavity 
in the presence of magnetic field perpendicular to the lid surface and its direction 
of motion; Re = 100 and Ha = 45. 




°V 0.2 0.4 0.6 0.8 1 

x/W 

Fig. 20. Effect of magnetic Prandtl number Pr^ on computed velocity profile u 
along the line y = W/2 and z = W/2 for simulation of 3D cubical lid-driven cavity 
in the presence of magnetic field perpendicular to the lid surface and its direction 
of motion; Re = 100 and Ha = 45. 
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Fig. 21. Computed velocity field u along the plane z = W/2 showing side wall-jets in 
a 3D cubical lid-driven cavity in the presence of magnetic field perpendicular to the 
lid surface and its direction of motion; Re = 100, Ha = 45 and Pr^ = 5.6 x 10"''. 

6 Summary and Conclusions 

In this paper, we devised a preconditioning approach for accelerating the 
steady state convergence of the solution of the generalized lattice Boltzmann 
equation (GLBE) with forcing term representing non-uniform force fields. A 
preconditioning parameter is introduced in the equilibrium moments and in 
the forcing terms of the GLBE that alleviates the disparity between the prop- 
agation speeds of density perturbation and fluid flow processes at low Mach 
numbers. The use of multiple relaxation times in the collision step involving 
the solution of the GLBE significantly improves the numerical stability of the 
approach as compared with the single relaxation time approach. Particular 
attention is paid to consistently preconditioning the projections of the forcing 
terms in the natural moment space of the GLBE. In particular, it is shown that 
to avoid spurious effects due to preconditioning, the slower processes involving 
the interaction of the velocity field and the external forces need to be precon- 
ditioned by the square of the prconditioning parameter. The limiting form of 
the consistent forcing term, which avoids discrete lattice artifacts particularly 
for spatially and temporally dependent external forces, is also obtained in the 
case of SRT-LBE. 

As an example of preconditioning an extended system of LBE for complex 
flows, we also developed a preconditioning procedure for a vector lattice ki- 
netic scheme that accelerate convergence of magnetohydrodynamics (MHD) 
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equations. In the case of MHD flows, there is an additional disparity between 
the propagation speeds of the perturbation of magnetic induction and the 
perturbation of the density field, which is mitigated by preconditioning an ad- 
ditional parameter. Motivated by applications in fusion engineering involving 
handling of liquid metals as blanket walls, in the context of preconditioning, 
we also devised a simple strategy to simulate flows with low magnetic Reynolds 
numbers or low magnetic Prandtl numbers. The greatest reduction in the num- 
ber of time steps to reach steady state is obtained when the preconditioning 
parameters of both the GLBE and the lattice kinetic scheme are the same. 
The preconditioned approach yielded solutions in very good agreement with 
prior solutions for several canonical examples. It is found that the precondi- 
tioned LBE reduced the number of steps for steady state convergence between 
several factors and as much as one or two orders of magnitude depending 
on the problem. In particular, for 3D MHD cavity flows that are character- 
ized by shear and other complex flow features, the preconditioning approach 
also yields accurate results in comparison with other numerical solutions. The 
preconditioning approach presented in this paper can be readily extended for 
other complex problems such as multiphase flows. 
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A Moments, Equilibrium Moments, Moment-projections of Source 
Terms 

A.l D2Q9 Model 

The components of the various elements in the moments are as follows [36]: 
/o = p, /i = e, /2 = e^, /s = ja;, /4 = qx, h = jy, /e = fr = Pj;^;, /s = Pxy 
Here, p is the density, e and represent kinetic energy that is independent of 
density and square of energy, respectively; and jy are the components of the 
momentum, i.e. = pu^ and jy = pUy, q^, and qy are the components of the 
energy flux, and p^x and Pxy are the components of the symmetric traceless 
viscous stress tensor. 

The corresponding components of the equilibrium moments, which are func- 
tions of the conserved moments, i.e. density p and momentum j , are as 
follows [36]: 

/^^ = p, ^ = -2p + 32i7 , ^ e^.e. ^ ^ _ ^7^^ fe, ^ fe, ^ 

The components of the source terms in moment space are functions of external 
force and velocity fields it , and are given as follows: 5*0 = 0, 5*1 = &{FxUx + 
FyUy),S2 — —6{FxUx + FyUy),Ss — Fx , S 4 — — Fx , S ^ — FyjSe — —Fy,S-j — 

2{FxUx - FyUy), Sg = {FxUy + FyUx)- 

A. 2 D3Q19 Model 

The components of the various elements in the moments are as follows [37]: 

/o = P, /i = e, /2 = /s = jx, fi = qx,h= jy, h = Qy, h = Jj, fs = Qzjg = 

'^Pxx, flO = '^'^xx, fu = Pww, f 12 = ^ww, fl3 = Pxy, f 14 = Pyz, fl5 = Pxz, f 16 ~ 

''^x, fi7 — "iTT-y, fi8 — "^z- Here, p is the density, e and represent kinetic 
energy that is independent of density and square of energy, respectively; jx, 

jy and jz are the components of the momentum, i.e. jx = pUx, jy = pUy, 
jz = pUz, Qx, Qy, Qz ai'c the components of the energy flux, and pxx, Pxy, Pyz 
and Pxz are the components of the symmetric traceless viscous stress tensor; 
The other two normal components of the viscous stress tensor, pyy and Pzz, 
can be constructed from Pxx and Pww, where p^jy, = Pyy —pzz- Other moments 
include -Rxx, t^ww, i^x, i^y and m^. The first two of these moments have the 
same symmetry as the diagonal part of the traceless viscous tensor pij, while 
the last three vectors are parts of a third rank tensor, with the symmetry of 

JkPmn- 
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The corresponding components of the equiUbrium moments, which are func- 
tions of the conserved moments, i.e. density p and momentum j , are as 
follows [37]: 



3Ja;) J< 



<i7 



3p; 



0, /rr^ = 



3iJ 
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T-7 



l„eg teq 
" ) J 13 
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0, 



a f=9 — „eg _ _2 ■ fcq . 
— -f^xj/ p ' J 14 — yyz 



--3p 



2 p 
'^1 = ^eg 



f eg _ ■ xeg _ 
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— 3Jz,J9 — 



peg 
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eg 
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peg 



/■eg _ 
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_ jxjz feq 
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The components of the source terms in moment space are functions of external 
force ^ and velocity fields Tt, and are given as follows [46]: 

5*0 = 0, S*i = SSiF^Ux + F^Uy + F^u^),S2 = -ll^F^u^ + FyUy + F^u^), S3 = 
Fx, 84^ = —^Fx, S^ = Fy, Sq = —^Fy, Sj = Fz, Ss = —^F^/Sg = 2{2FxUx — 
FyUy - F^Uz),Sio^= -{2FxUx - FyUy -J^zUz),Sn = 2{FyUy - F^Uz),Si2 = 

-{FyUy - F;,Uz),Sl3 = {F^Uy + FyUx).,Su = {FyU z + F;,Uy),Sl5 = (F^Uz + 

F,u,),Sie^O,Sn^O,Sis^O. 



B Chapman-Enskog Analysis of the Preconditioned GLBE with 
Forcing Term for D2Q9 Model 



In this section, by employing the Chapman-Enskog multiscale analysis [48,39], 
we derive the macroscopic dynamical equations for the preconditioned GLBE 
with forcing term corresponding to the D2Q9 model. The analysis of the pre- 
conditioned GLBE with other two- or three-dimensional models can be carried 
out in an analogous way. First, we introduce the expansions 

00 

f=^e"f("), (B.l) 

n=0 
00 

9*=E^"^*n, (B.2) 

n=0 

where e — 5t, along with the Taylor series into the preconditioned GLBE pre- 
sented in Sec. 3. Then, recognizing that it was derived after making use of the 

the transformation f = f — l/2Sona second-order time discretization of the 
source terms to make it effectively explicit, and dropping the "overbar" sub- 
sequently for convenience, the following equations are obtained as consecutive 
orders of the parameter e in moment space (see Ref. [39] for details) 
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0(e°) :f(°) (B.3) 
O(e^) : [dt, + m) f = -A*f(i) + S*, (B.4) 

o{e') : + + (i - h*^ f^^) = -A*f^'\ (b.s) 

where = TeaiT~^. As it is known that the non-hnear terms in the equihb- 
rium moments give rise to slower convective behavior of the fluid as compared 
with the propagation speed of the density perturbations, we precondition those 
terms by the parameter 7 to obtain f"^"?'*. The forcing terms in moment space 
are functions of external forces and velocity fields as given in Appendix A. We 
assume that for the components linear in T^, we precondition them by 7. On 
the other hand, we precondition those that are non-linear due to interactions 
between ^ and it by an unknown parameter 71, whose form will be deduced 
as part of the analysis. The final expressions for the preconditioned quantities 
f^^'* and S*, as well as A* are given in Sec. 3. 

The components of the flrst-order equations in moment space, i.e. Eq. (B.4 
can be written as 

dtoP + dxjx + dyjy = (B.6 
dt, (-2p + -^"j = -ste(i) + (B.7 

V 7P y 71 

dt, (p - 3^-^] - ^ • 7 = -«;e^W - 6^ (B.S 
\ IP ) 7i 

dt^j. + i\p +^]+dy f ^) = ^ (B.9 

V> IP) \1P ) 7 

a. (-i) + 4 (4p - ^^f) + a, (^) = - ^ (B.10 
0^ i-j,) + (^) + a, (-ip + ^--i ) = - Si (B.12 
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Similarly, the components of the second-order equations in moment space, i.e. 
Eq. (B.5) are 

dt,p = (B.15) 



-2p + 



1 - ^s: 



IP 



'to 



= (1) 



+ 



1 - -st 



.(1) 



(2) 



(B.16) 
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(B.17) 



dtjx + ( g 
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ft, 



1 - -si 
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(B.18) 



^1 (-Jx) + di 
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(B.19) 



dtjy + ft 
1 

6 



1 - -si 

2 ^ 



AD 

fxy 



ft 



=(1) _ 



(B.20) 



fti (-Jy) + ft 



ft, 
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1 st 

2 ^. 



e(i) + - 



1 - 



dr 



l-^sl 



fxy 



+ 



e2(i) + 



1 - 



-«6??lB.21) 
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(B.22) 
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Jxjy 
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1- 



(B.23) 



To obtain preconditioned hydrodynamical equations, we need to combine the 
evolution equations of moments corresponding to Eqs. (B.6), (B.9) and (B.ll) 
with Eqs. (B.15), (B.18) and (B.20), respectively. Inspection of these six equa- 
tions show that we need explicit expressions for the following non-equilibrium 
moments: e^^\ andpg. Thus, from Eqs. (B.7), (B.13) and (B.14), we get 



.(1) = ± 
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-Si 



to 



-2p + 3 
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IP 



+ 6- 



u 



71 



fxx 



yxy 
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'fx 



to 



to 



- dx l^Jx ) +dy{-Jy) + 2 



(B.24) 

Fri-Ux FyUy) 



IP 



d:r. 



a, 



Jy + 



7i 
7i 



B.25) 



(B.26) 



Eqs. (B.24)-(B.26) require time derivatives of the density and velocity (or 
momentum) fields, which can be obtained from Eqs. (B.6), (B.9) and (B.ll) 
and truncating terms of O(Ma^) or higher. Thus, we have 



iFxUx 



*° V P / ~ 



(FxUy+FyUx) 
7 



Substituting the above relations in Eq. (B.24) 



3 I _^ 2_ 

7 



u 



7i 



(B.27) 



In order to eliminate the dependence of forcing terms in the above equation, 
Eq. (B.27), so that the macrodynamical equations recover correct physics with- 
out any spurious effects, we need to set 

71 = 7'- (B.28) 
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Thus, Eq. (B.28) suggests that the moment projections of forcing terms in- 
volving non-hnear interactions of external force and velocity fields should be 
preconditioned by a factor of 1/7^, as they represent slower physical processes 
than fluid flow itself. Hence, we get 

e« = -2!^ . t. (B.29) 



Similarly, 



P^^I - -~ (d^Jx - dyjy) , (B.30) 



and 



1 1 



{dxjy + dyj:c) 



(B.31) 



So, the preconditioned dynamical equations of the conserved moments are fi- 
nally obtained by adding Eqs. (B.6), (B.9) and (B.ll) to Eqs. (B.15), (B.18) 
and (B.20), respectively, after multiplying the latter with 5t, and using Eqs. (B.29)- 
(B.31). They correspond to the following preconditioned weakly compressible 
Navier-Stokes equations 

dtP + dJ, + dyjy = Q (B.32) 
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where the pressure field p is given by 



(B.35) 



and the transport coefficients, viz., the bulk and shear viscosities, respectively, 
as 



1 1 



(B.36) 



and 



1/ — 'f- 



4 2 



St, /3 = 7,i 



(B.37) 



C Chapman-Enskog Analysis of the Preconditioned Vector Ki- 
netic Equation 



We now perform the Chapman-Enskog analysis of the preconditioned vector 
kinetic equation by introducing the expansions 



a ^Ve^o^"*^ 

n=0 
oo 



(C.l) 
(C.2) 



n=0 



with e = St, in conjunction with the Taylor series [48,51], which result in the 
following: 



(0) 
ai 



^ (1) 

m 



(C.3) 
(C.4) 

(C.5) 



Now, using the following summational constraints I]^=o 9aj — ^jt Sa=o ^ai9aj 
Ag\ T!'-=o9^:^ = and E^=o e..^?if = A^f for n > 1, with = 
and taking zeroth moments of Eqs. (C.4) and (C.5), we get 
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2t* 



0. 



(C.6) 
(C.7) 



Taking the first moment, i.e. J2^c^o^aj{-) of Eq. (C.4) and using the identity 



Ea=o ^ajGakQaP = ^m^jkBi, where 5jk is the Kronecker delta, we get 



3% m 



(C.8) 



With the scaling 0{Bi) ~ 0{ui) [51] and considering the zeroth-order momen- 
tum and magnetic induction equations, dtgAj^^ ~ O(Ma^) and hence can be 
neglected. As a result, 

^ -r*JmdiB,. (C.9) 



Finally, adding Eq. (C.6) and Eq. {Cl)x5t and using dt ~ dt^ +5tdti, along 
with Eq. (C.9), we get the preconditioned magnetic induction equation 

dtBi + —Vj {ujBi - BjUi) = —dj ivdjBi) , (C. 10) 



where 77 = -fmOm {r^ - 1/2) 5*. 
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